Trong model sẽ bao gồm 2 phần cần chú ý:
Quan trọng nhất là model “tốt nhất” (theo tuỳ chuẩn) không có nghĩa là model đúng với thực tế “true”
All models are wrong, but some are useful. - George Box
For such a model there is no need to ask the question “Is the model true?”. If “truth” is to be the “whole truth” the answer must be “No”. The only question of interest is “Is the model illuminating and useful?”.
Đầu tiên load các libary cần thiết
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.1 ✔ purrr 1.0.1
## ✔ tibble 3.1.8 ✔ dplyr 1.1.0
## ✔ tidyr 1.3.0 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 1.0.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(modelr)
options(na.action = na.warn)
Đầu tiên xem data set của sim1 trong model r, gồm 2 biến x và y ở dạng continuous
sim1
## # A tibble: 30 × 2
## x y
## <int> <dbl>
## 1 1 4.20
## 2 1 7.51
## 3 1 2.13
## 4 2 8.99
## 5 2 10.2
## 6 2 11.3
## 7 3 7.36
## 8 3 10.5
## 9 3 10.5
## 10 4 12.4
## # … with 20 more rows
ggplot(sim1, aes(x, y)) +
geom_point()
Từ plot có thể thấy hai biến x và y có correlation mạnh và ở dạng tuyến tính, y = a_0 + a_1*x.
Thử dùng geom_abline để vẽ các đường thẳng với intercept và slope bất kỳ
models <- tibble(
a1 = runif(250, -20, 40),
a2 = runif(250, -5, 5)
)
ggplot(sim1, aes(x, y)) +
geom_abline(aes(intercept = a1, slope = a2), data = models, alpha = 1/4) +
geom_point()
=> Phải tìm giá trị a_0 và a_1 khớp nhất với data.
Sử dụng khoảng cách khác biệt giá trị y dự đoán và giá y trị thực tế để tìm model thích hợp nhất.
Ví dụ, tính giá trị y dự đoán nếu a_0 = 7 và a_1 =1.5
y = a_0 + a_1 *x
y = 7 + 1.5 *x
model1 <- function(a, data) {
a[1] + data$x * a[2]
}
model1(c(7, 1.5), sim1) #y = 7 + 1.5 *x
## [1] 8.5 8.5 8.5 10.0 10.0 10.0 11.5 11.5 11.5 13.0 13.0 13.0 14.5 14.5 14.5
## [16] 16.0 16.0 16.0 17.5 17.5 17.5 19.0 19.0 19.0 20.5 20.5 20.5 22.0 22.0 22.0
viết function để tính độ lệch
measure_distance <- function(mod, data) {
diff <- data$y - model1(mod, data)
sqrt(mean(diff ^ 2))
}
measure_distance(c(7, 1.5), sim1)
## [1] 2.665212
Dùng purrr package để tính căn bậc hai của sai số bình phương cho 250 models
sim1_dist <- function(a1, a2) {
measure_distance(c(a1, a2), sim1)
}
models <- models %>%
mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
models
## # A tibble: 250 × 3
## a1 a2 dist
## <dbl> <dbl> <dbl>
## 1 -6.62 3.94 5.86
## 2 22.2 0.186 9.60
## 3 11.9 -0.0713 7.62
## 4 -0.692 2.61 3.23
## 5 14.8 -1.29 12.6
## 6 32.0 1.86 26.9
## 7 6.70 3.95 14.2
## 8 -19.9 -1.88 47.2
## 9 2.43 -2.86 32.1
## 10 -12.6 -0.293 30.6
## # … with 240 more rows
Thử vẽ 10 model với độ lệch nhỏ nhất
ggplot(sim1, aes(x, y)) +
geom_point(size = 2, colour = "grey30") +
geom_abline(
aes(intercept = a1, slope = a2, colour = -dist),
data = filter(models, rank(dist) <= 10)
)
Ngoài ra còn có thể chuyển giá trị a_1 a_2 vào scatterplot để xem
ggplot(models, aes(a1, a2)) +
geom_point(data = filter(models, rank(dist) <= 10), size = 4, colour = "red") +
geom_point(aes(colour = -dist))
Bỏ scatterplot này vào lưới để chọn giá trị a1 và a2 một cách đồng đều hơn
grid <- expand.grid(
a1 = seq(-5, 20, length = 25),
a2 = seq(1, 3, length = 25)
) %>%
mutate(dist = purrr::map2_dbl(a1, a2, sim1_dist))
grid %>%
ggplot(aes(a1, a2)) +
geom_point(data = filter(grid, rank(dist) <= 10), size = 4, colour = "red") +
geom_point(aes(colour = -dist))
Đem 10 models theo scatter plot này vào lại data
ggplot(sim1, aes(x, y)) +
geom_point(size = 2, colour = "grey30") +
geom_abline(
aes(intercept = a1, slope = a2, colour = -dist),
data = filter(grid, rank(dist) <= 1)
)
Nếu làm cho grid càng nhỏ thì cuối cùng ta sẽ có model fit “best”.
Hoặc sử dụng phương pháp Newton-Raphson: chọn một điểm rồi tìm slope dốc nhất.
best <- optim(c(0, 0), measure_distance, data = sim1)
best$par
## [1] 4.222248 2.051204
ggplot(sim1, aes(x, y)) +
geom_point(size = 2, colour = "grey30") +
geom_abline(intercept = best$par[1], slope = best$par[2])
Phương pháp này ứng dụng được trên mọi dạng model, nhưng trong model linear, có function lm để fit model tốt nhất với data
sim1_mod <- lm(y ~ x, data = sim1)
coef(sim1_mod)
## (Intercept) x
## 4.220822 2.051533
Tạo data dự đoán
grid <- sim1 %>%
data_grid(x)
grid
## # A tibble: 10 × 1
## x
## <int>
## 1 1
## 2 2
## 3 3
## 4 4
## 5 5
## 6 6
## 7 7
## 8 8
## 9 9
## 10 10
Thêm giá trị y dự đoán vào
sim1
## # A tibble: 30 × 2
## x y
## <int> <dbl>
## 1 1 4.20
## 2 1 7.51
## 3 1 2.13
## 4 2 8.99
## 5 2 10.2
## 6 2 11.3
## 7 3 7.36
## 8 3 10.5
## 9 3 10.5
## 10 4 12.4
## # … with 20 more rows
#y = a_1 + a_2 * x
grid <- grid %>%
add_predictions(sim1_mod)
grid
## # A tibble: 10 × 2
## x pred
## <int> <dbl>
## 1 1 6.27
## 2 2 8.32
## 3 3 10.4
## 4 4 12.4
## 5 5 14.5
## 6 6 16.5
## 7 7 18.6
## 8 8 20.6
## 9 9 22.7
## 10 10 24.7
Vẽ data dự đoán vào trong plot data quan sát
ggplot(sim1, aes(x)) +
geom_point(aes(y = y)) +
geom_line(aes(y = pred), data = grid, colour = "red", linewidth = 1)
Function add_residuals() để thêm giá trị residual vào bảng prediction
sim1 <- sim1 %>%
add_residuals(sim1_mod)
sim1
## # A tibble: 30 × 3
## x y resid
## <int> <dbl> <dbl>
## 1 1 4.20 -2.07
## 2 1 7.51 1.24
## 3 1 2.13 -4.15
## 4 2 8.99 0.665
## 5 2 10.2 1.92
## 6 2 11.3 2.97
## 7 3 7.36 -3.02
## 8 3 10.5 0.130
## 9 3 10.5 0.136
## 10 4 12.4 0.00763
## # … with 20 more rows
Cách phương pháp vẽ graph cho residual
ggplot(sim1, aes(resid)) +
geom_freqpoly(binwidth = 0.5)
Hoặc ở dạng scatterplot
ggplot(sim1, aes(x, resid)) +
geom_ref_line(h = 0) +
geom_point()
Scatterplot của resid ở dạng random noise -> không có bias
plot(sim1_mod)
Dùng model_matrix() để xem dạng model
df <- tribble(
~y, ~x1, ~x2,
4, 2, 5,
5, 1, 6
)
df
## # A tibble: 2 × 3
## y x1 x2
## <dbl> <dbl> <dbl>
## 1 4 2 5
## 2 5 1 6
model_matrix(df, y ~ x1)
## # A tibble: 2 × 2
## `(Intercept)` x1
## <dbl> <dbl>
## 1 1 2
## 2 1 1
dạng model là: y ~ x1 tương đương với y =
Thêm -1 để bỏ cột intercept mặc định
model_matrix(df, y ~ x1 - 1)
## # A tibble: 2 × 1
## x1
## <dbl>
## 1 2
## 2 1
y = 2 + x
model_matrix(df, y ~ x1 + x2)
## # A tibble: 2 × 3
## `(Intercept)` x1 x2
## <dbl> <dbl> <dbl>
## 1 1 2 5
## 2 1 1 6
Model cơ bản y ~ x viết ở dưới dạng y = a_1 + a_2 * x
Để xem cách R làm model, dùng function model_matrix
df <- tribble(
~y, ~x1, ~x2,
4, 2, 5,
5, 1, 6
)
model_matrix(df, y ~ x1)
## # A tibble: 2 × 2
## `(Intercept)` x1
## <dbl> <dbl>
## 1 1 2
## 2 1 1
R tự động thêm cột intercept vào model với cột đầu chỉ có số 1. Cách để bỏ cột số một trong model_matrix
model_matrix(df, y ~ x1 -1 )
## # A tibble: 2 × 1
## x1
## <dbl>
## 1 2
## 2 1
model_matrix(df, y ~ x1 + x2)
## # A tibble: 2 × 3
## `(Intercept)` x1 x2
## <dbl> <dbl> <dbl>
## 1 1 2 5
## 2 1 1 6
Tạo data frame mẫu và model response phụ thuộc vào giới tính
df <- tribble(
~ sex, ~ response,
"male", 1,
"female", 2,
"male", 1
)
model_matrix(df, response ~ sex)
## # A tibble: 3 × 2
## `(Intercept)` sexmale
## <dbl> <dbl>
## 1 1 1
## 2 1 0
## 3 1 1
sim2 %>% head()
## # A tibble: 6 × 2
## x y
## <chr> <dbl>
## 1 a 1.94
## 2 a 1.18
## 3 a 1.24
## 4 a 2.62
## 5 a 1.11
## 6 a 0.866
ở dạng graph
ggplot(sim2) +
geom_point(aes(x, y))
fit model vào
mod2 <- lm(y ~ x, data = sim2)
grid <- sim2 %>%
data_grid(x) %>%
add_predictions(mod2)
grid
## # A tibble: 4 × 2
## x pred
## <chr> <dbl>
## 1 a 1.15
## 2 b 8.12
## 3 c 6.13
## 4 d 1.91
Prediction của y ở mỗi x chính là trung bình của y ở điểm đó
ggplot(sim2, aes(x)) +
geom_point(aes(y = y)) +
geom_point(data = grid, aes(y = pred), colour = "red", size = 4)
Trong trường hợp sim 3 vừa có continuous and categorical
ggplot(sim3, aes(x1, y)) +
geom_point(aes(colour = x2))
2 models có thể fit cho data này
mod1 <- lm(y ~ x1 + x2, data = sim3)
mod2 <- lm(y ~ x1 * x2, data = sim3)
(y ~ x1 + x2): y = a_0 + a_1*x1 + a_2*x2
(y ~ x1 * x2): y = a_0 + a_1*x1 + a_2*x2 +a_12*x1*x2
Tạo bảng prediction
grid <- sim3 %>%
data_grid(x1, x2) %>%
gather_predictions(mod1, mod2)
grid
## # A tibble: 80 × 4
## model x1 x2 pred
## <chr> <int> <fct> <dbl>
## 1 mod1 1 a 1.67
## 2 mod1 1 b 4.56
## 3 mod1 1 c 6.48
## 4 mod1 1 d 4.03
## 5 mod1 2 a 1.48
## 6 mod1 2 b 4.37
## 7 mod1 2 c 6.28
## 8 mod1 2 d 3.84
## 9 mod1 3 a 1.28
## 10 mod1 3 b 4.17
## # … with 70 more rows
Graph
ggplot(sim3, aes(x1, y, colour = x2)) +
geom_point() +
geom_line(data = grid, aes(y = pred)) +
facet_wrap(~ model)
Xem residual của từng line
sim3 <- sim3 %>%
gather_residuals(mod1, mod2)
ggplot(sim3, aes(x1, resid, colour = x2)) +
geom_point() +
facet_grid(model ~ x2)
mod1 <- lm(y ~ x1 + x2, data = sim4)
mod2 <- lm(y ~ x1 * x2, data = sim4)
grid <- sim4 %>%
data_grid(
x1 = seq_range(x1, 5),
x2 = seq_range(x2, 5)
) %>%
gather_predictions(mod1, mod2)
grid
## # A tibble: 50 × 4
## model x1 x2 pred
## <chr> <dbl> <dbl> <dbl>
## 1 mod1 -1 -1 0.996
## 2 mod1 -1 -0.5 -0.395
## 3 mod1 -1 0 -1.79
## 4 mod1 -1 0.5 -3.18
## 5 mod1 -1 1 -4.57
## 6 mod1 -0.5 -1 1.91
## 7 mod1 -0.5 -0.5 0.516
## 8 mod1 -0.5 0 -0.875
## 9 mod1 -0.5 0.5 -2.27
## 10 mod1 -0.5 1 -3.66
## # … with 40 more rows
Sử dụng seq_range() thay thế cho data_grid(). Thay thế cho việc sử dụng giá trị unique của x
Sử dụng các argument như
Ví dụ
seq_range(c(0.0123, 0.923423), n = 5)
## [1] 0.0123000 0.2400808 0.4678615 0.6956423 0.9234230pretty = TRUE
seq_range(c(0.0123, 0.923423), n = 5, pretty = TRUE)
## [1] 0.0 0.2 0.4 0.6 0.8 1.0trim = 0.1, cắt 10% giá trị đuôi
x1 <- rcauchy(100)
print(x1)
## [1] 7.17455599 0.51857630 1.50053731 -173.50882444 10.55462922
## [6] -0.12366276 0.03673073 -2.05993929 3.14016643 2.87335478
## [11] 1.17228076 -1.78004277 -1.64912075 0.37686329 -0.93382796
## [16] -1.41491860 -1.27807239 0.60098207 1.15952963 -1.37894904
## [21] 0.34351244 -15.64513911 -5.17460388 -0.67602844 4.66956926
## [26] -0.06596813 -0.42246329 0.66521295 0.45269676 0.46496412
## [31] 4.80513029 -2.58605494 -5.32788465 -3.30938897 -1.68157543
## [36] 3.77587123 0.58505480 0.30220331 0.43743750 -0.01776576
## [41] -0.09188035 10.58232172 -2.87903719 -0.84014695 8.91803085
## [46] -2.32556547 -4.61932400 -0.25750880 -0.31349105 0.29792883
## [51] 2.83804175 2.08793754 15.70480489 1.30658619 -3.13945941
## [56] 2.97463983 -1.68151356 -1.16133722 -0.71174224 -0.71249720
## [61] 0.41676315 0.20540708 0.50853307 -0.62810086 -1.16936187
## [66] 1.88497384 2.08855329 -0.96490480 0.15635953 2.08706652
## [71] 1.62181829 -0.34904277 1.48756045 -0.19601606 -0.56082780
## [76] -0.18970034 -0.93878887 -25.24477495 8.03804283 5.57527568
## [81] -5.46078151 -0.41090739 -0.12711524 2.37673550 8.02653473
## [86] -3.11376245 -1.01051437 4.87194781 -1.47955698 0.41531512
## [91] 1.29209637 -1.90283696 -6.13795523 0.25426971 -6.46131391
## [96] -1.28506105 -0.81624626 -0.03923384 0.64624844 2.45024668
seq_range(x1, n = 5, trim = 0.1)
## [1] -5.494640 -2.114203 1.266235 4.646673 8.027110vẽ
ggplot(grid, aes(x1, x2)) +
geom_tile(aes(fill = pred)) +
facet_wrap(~ model)
Không thấy rõ model nào fit tốt hơn
ggplot(grid, aes(x1, pred, colour = x2, group = x2)) +
geom_line() +
facet_wrap(~ model)
ggplot(grid, aes(x2, pred, colour = x1, group = x1)) +
geom_line() +
facet_wrap(~ model)
log(y) ~ sqrt(x1) + x2 biến đổi thành log(y) = a_1 + a_2 * sqrt(x1) + a_3 * x2
Nếu transformation sử dụng dấu +, *, ^ hay -, ta phải bỏ vào I( )
y ~ x + I(x ^ 2) biến đổi thành y = a_1 + a_2 * x + a_3 * x^2
R sẽ tự động bỏ biến dư nếu viết dạng x + x
y ~ x ^ 2 + x sẽ
thành y = a_1 + a_2 * x
df <- tribble(
~y, ~x,
1, 1,
2, 2,
3, 3
)
model_matrix(df, y ~ x^2 + x)
## # A tibble: 3 × 2
## `(Intercept)` x
## <dbl> <dbl>
## 1 1 1
## 2 1 2
## 3 1 3
model_matrix(df, y ~ I(x^2) + x)
## # A tibble: 3 × 3
## `(Intercept)` `I(x^2)` x
## <dbl> <dbl> <dbl>
## 1 1 1 1
## 2 1 4 2
## 3 1 9 3
Chọn polynomial để smooth tốt
model_matrix(df, y ~ poly(x, 2))
## # A tibble: 3 × 3
## `(Intercept)` `poly(x, 2)1` `poly(x, 2)2`
## <dbl> <dbl> <dbl>
## 1 1 -7.07e- 1 0.408
## 2 1 -9.07e-17 -0.816
## 3 1 7.07e- 1 0.408
Khi sử dụng poly, ngoài range của data, polynomial sẽ tự động dự đoán positive and negative infinity. Có sử dụng natural spline
library(splines)
model_matrix(df, y ~ ns(x, 2))
## # A tibble: 3 × 3
## `(Intercept)` `ns(x, 2)1` `ns(x, 2)2`
## <dbl> <dbl> <dbl>
## 1 1 0 0
## 2 1 0.566 -0.211
## 3 1 0.344 0.771
sim5 <- tibble(
x = seq(0, 3.5 * pi, length = 50),
y = 4 * sin(x) + rnorm(length(x))
)
ggplot(sim5, aes(x, y)) +
geom_point()
Fit thử 5 loại spline model
mod1 <- lm(y ~ ns(x, 1), data = sim5)
mod2 <- lm(y ~ ns(x, 2), data = sim5)
mod3 <- lm(y ~ ns(x, 3), data = sim5)
mod4 <- lm(y ~ ns(x, 4), data = sim5)
mod5 <- lm(y ~ ns(x, 5), data = sim5)
grid <- sim5 %>%
data_grid(x = seq_range(x, n = 50, expand = 0.1)) %>%
gather_predictions(mod1, mod2, mod3, mod4, mod5, .pred = "y")
ggplot(sim5, aes(x, y)) +
geom_point() +
geom_line(data = grid, colour = "red") +
facet_wrap(~ model)
df <- tribble(
~x, ~y,
1, 2.2,
2, NA,
3, 3.5,
4, 8.3,
NA, 10
)
mod <- lm(y ~ x, data = df)
## Warning: Dropping 2 rows with missing values
mod <- lm(y ~ x, data = df, na.action = na.exclude)
nobs(mod)
## [1] 3